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The non-equilibrium relaxational properties of a three dimensional Coulomb glass model are in- 
vestigated by kinetic Monte Carlo simulations. Our results suggest a transition from stationary 
to non- stationary dynamics at the equilibrium glass transition temperature of the system. Below 
the transition the dynamic correlation functions loose time translation invariance and electron dif- 
fusion is anomalous. Two groups of carriers can be identified at each time scale, electrons whose 
motion is diffusive within a selected time window and electrons that during the same time interval 
remain confined in small regions in space. During the relaxation that follows a temperature quench 
an exchange of electrons between these two groups takes place and the non-equilibrium excess of 
diffusive electrons initially present decreases logarithmically with time as the system relaxes. This 
bimodal dynamical heterogeneity persists at higher temperatures when time translation invariance 
is restored and electron diffusion is normal. The occupancy of the two dynamical modes is then 
stationary and its temperature dependence reflects a crossover between a low-temperature regime 
with a high concentration of electrons forming fluctuating dipoles and a high-temperature regime 
in which the concentration of diffusive electrons is high. 
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I. INTRODUCTION 

Recent experimental studies of hopping conductance 
in Anderson insulators showed striking non-equilibrium 
effects that persist for long times at low temperature. 
i&S&k These results support early theoretical predic- 
tions of a low-temperature glassy phase in interacting dis- 
ordered electron systems in the strongly localized limit. 6,7 
In this regime the essential physics is well captured by the 
Coulomb glass model that describes a system of interact- 
ing electrons hopping between randomly distributed sites 
that correspond to the localization centers of the single- 
electron wavefunctions&LSiSiiSiii 

Although Coulomb glass models have been extensively 
studied during the last thirty years and their equilibrium 
properties are by now fairly well understood&2i2i2ii24iiii^ 
much less is known about the properties of Coulomb 
glasses out of equilibriumiiiiiiiSiiii 

In Reference ^| the off-equilibrium dynamics of the 
Coulomb glass was investigated by studying the scaling 
properties of the non-stationary correlation and response 
functions, a tool often used in the study of other glassy 
systems such as structural and spin glasses jiLigiiiLSii . It 
has recently been realized that further insights on the 
nature of the glassy phase can be gained from the analy- 
sis of dynamical heterogeneitiesiSi* 2 ^ 2 ^ This approach is 
particularly appealing in the context of Coulomb glasses 
since the low-temperature hopping conductance is dom- 
inated by the diffusion of carriers on a set of conducting 



percolation paths&i We expect interactions between the 
electrons to have strong effects on this type of motion 
since hops of individual carriers alter the effective ran- 
dom potential felt by the others and thus modify the 
structure of the percolating network. Fluctuation effects 
of this type were shown to lead to low-frequency l//-likc 
noise in the conductance of Coulomb glasses. 24,25 

In this paper we investigate the dynamical proper- 
ties of the three dimensional random-site Coulomb glass 
model by kinetic Monte Carlo simulation using a realistic 
microscopic dynamics that favors the emergence of local 
effective constraints in the kinetics. In our simulations 
the system is initially quenched from infinite temperature 
to a working temperature T and its evolution in time is 
characterized for different values of T through the time 
dependence of the relevant correlation functions. 

One of our main observations is the appearance of a 
dynamical crossover from equilibrium dynamics to slow 
non-equilibrium dynamics at a temperature T g ~ T c , 
where T c is the equilibrium freezing transition temper- 
ature of the modeliii This crossover takes place even for 
relatively small system sizes and occurs at the temper- 
ature at which the equilibration time of the finite sam- 
ple becomes much longer than the time-scale of the sim- 
ulation. In this regime the time-dependent correlation 
functions exhibit slow relaxation and have aging proper- 
ties. We found that the dynamics of the Coulomb glass 
is heterogeneous as observed in other glassy systems 
Heterogeneities can be characterized by examination of 
the evolution of diffusion fronts. A statistical analysis 
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of the electron trajectories over a fixed time window 
shows that most carriers belong to one or the other of 
two groups. The first group is that of those electrons 
that diffused away from their initial position during the 
chosen time interval. The second group is that of the 
electrons that remained confined in relatively small re- 
gions in space during that time. An exchange of electrons 
between these two dynamical modes takes place as the 
system relaxes after a temperature quench. In the aging 
regime this exchange is very slow and the non-equilibrium 
excess of electrons with metallic hopping present in the 
system right after the quench decreases logarithmically 
with time. In this temperature range the diffusion is 
anomalous. We found that the mean squared displace- 
ment (Ax 2 (t)) ~ tf 1 where the exponent is rj < 1 and 
depends on the age of the system. In the equilibrium 
regime time translation invariance is recovered and we 
observe normal diffusion. However, the bimodal dynam- 
ical heterogeneity still persists. The occupancy of the 
two dynamical modes is stationary and its temperature 
dependence reflects a crossover from a low-temperature 
regime with a high concentration of electrons forming 
fluctuating dipoles to a high-temperature regime in which 
the concentration of diffusive electrons is high. 

The paper is organized as follows. SectionlTllcontams a 
description of the model and of our numerical method. In 
Section UTTl we discuss the properties of the local-density 
autocorrelation function in and out of equilibrium. Sec- 
tion Hy| is devoted to the analysis of electron diffusion in 
the system. In Section we study the statistical prop- 
erties of the diffusion fronts and show that they provide 
evidence for the existence of heterogeneous transport in 
the system. Finally, we summarize the conclusions of our 
study in Section IVT1 



II. THE MODEL 

The Hamiltonian of the classical three dimensional 
Coulomb glass is£ 
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where Ri denotes center of localization of a single- 
particle localized electronic state, ipi the energy of the 
state and e and k are the electron charge and the 
medium's dielectric constant, respectively. Strong on-site 
correlations limit the occupancy of the electronic states 
to rii = 0,1. Charge neutrality is assured by a uniform 
compensating positive charge density K = (1/N) J2i n i- 
The positions Ri of the localized states and their ener- 
gies ipi are both random variables. However, it is common 
practice to study two complementary simplified versions 
of the model. These are respectively referred to as the 
lattice model and the random site model in Reference 
TTEL In the lattice model the sites are assumed to lie on 
a regular cubic lattice and only the randomness in the 



energies tpi is taken into consideration. Conversely, in 
the random site modeliSiii only the positional disorder 
is taken into account and tpi = 0. It is standard prac- 
tice to concentrate on the particle-hole symmetric case 
K — 1/2 for which the analysis of the results is simpler. 

While it has been established that the three dimen- 
sional random-site model has an equilibrium glass tran- 
sition at low temperature^ it is not yet clear whether 
this is also the case for the lattice model ^'-'ffi Therefore, 
we shall only discuss the dynamics of the random-site 
model in this paper. 

To model the dynamics of the system we let it evolve 
through sequential single-electron hops from occupied 
sites a to empty sites b. The transition rate, that mimics 
phonon assisted processes, is 



T a ^ b = To l e- 2R ^ min[l,e- AB -»/ T 
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where To is a microscopic time, £ is the spatial exten- 
sion of the localized wavefunctions and R a b = |R a — R&l- 
AE a b, the total energy difference in the transition, is 
given by 
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Srii = iii — K. 

The first factor in Eq.(J2J) reflects the exponential de- 
cay of the electron-phonon matrix element between two 
electronic wavefunctions centered at positions R a and 
Rb. The second factor is the thermal part of the transi- 
tion probability. In Monte Carlo simulations performed 
by other authors the transition probability is taken in- 
dependent of the distance R a b [the first exponential in 
the equivalent of Eq.(2) is absent] Aii This type of non- 
local dynamics, convenient for rapid equilibration, may 
not be appropriate for the study of off-equilibrium relax- 
ation. With the local dynamics of Eq. [21 electron hops 
that decrease the energy are essentially restricted to a 
region whose linear size is the localization length £. This 
introduces dynamic constraints that contribute to make 
the relaxation out of an excited configuration slower. 

In our simulations we take the mean distance be- 
tween sites ao as the unit of length, the Coulomb energy 
Ec = e 2 /(k<2o) as the unit of energy and we choose for 
convenience £ = ao- 

We simulated systems of N = L 3 sites and M = N/2 
electrons for samples with L =6, 8 and 10 in the temper- 
ature range 0.01 < T < 0.1. The localization centers are 
distributed randomly and uniformly inside a computa- 
tional cubic box of side L and we take periodic boundary 
conditions in all directions. To simulate a quench from 
high temperature we start from a random electron con- 
figuration at t — and let the system freely evolve with 
the dynamics (O at the working temperature T. The 
elementary Monte Carlo move consists of an attempt to 
move an electron from a randomly chosen occupied site 
a to an empty site b. Once a is chosen, the destination 
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site b is chosen randomly using the probability distribu- 
tion of the hoppings, P (R a b) °c ex P ( — 2i? ah ). A cutoff at 
-Rah = L/2 is imposed by our using of periodic boundary 
conditions and restricts us to a range of temperatures 
for which hops at distances R a b ~ L can be neglected. 
The probability of acceptance of the move is the ther- 
mal factor in Eq. @. If the hop is accepted, the vector 
displacement of the hopping electron <5r is defined as the 
vector going from site a to the closest periodic image of 
site b. A Monte Carlo step (MCS) consists of N hopping 
attempts. Our runs were typically 2 x 10 6 MCS long. 
Physical quantities were monitored as a function of time 
for each sample and the results were averaged over be- 
tween 150 and 600 realizations of the disorder and initial 
conditions, depending on the system size and tempera- 
ture. 



III. THE LOCAL DENSITY 
AUTOCORRELATION FUNCTION 

The two-time charge autocorrelation function isi& 

C(t + t w ,t w ) = jj22(6rii(t + t w ) 6m{t w )) , (4) 

i 

where the brackets denote the average over configura- 
tional disorder, initial conditions, and thermal noise, and 
the waiting time t w is the time elapsed since the quench 
from infinite temperature. 

The function C(t + t w , t w ) is the overlap of the charge 
configurations at times t + t w and t w . If t w is larger 
than the equilibration time r cq the state of the system 
is time-translational invariant and the correlation func- 
tion depends only on the time difference t. Otherwise, C 
depends on both t and t w . 

We describe in the following results for a system of lin- 
ear size L = 10. Fig.^displays C(t+t w , t w ) vs. t for three 
waiting times, t w — 10 3 , 10 4 and 10 5 and two represen- 
tative cases, T = 0.07 [FigEta)] and T = 0.03 [FigEtb)]. 
We observe that in the first case C(t + t w ,t w ) w C(t), i.e. 
the system is time translational invariant. This means 
that equilibrium was reached within a time shorter than 
the shortest of the waiting times considered, r eq (T) < 
10 3 . C{t) is thus the equilibrium relaxation function. In 
the second case, however, time translation invariance is 
lost and C(t + t w , t w ) depends on both time arguments, 
a phenomenon known as aging. The system was thus 
unable to equilibrate within the time scale of the simu- 
lation. Note that for each value of t w the relaxation is 
very slow (roughly logarithmic) and becomes slower with 
increasing t w . 

We found that non-equilibrium relaxation appears be- 
low a dynamic crossover temperature T g ~ 0.05. The 
value of T g is remarkably close to the equilibrium tran- 
sition temperature of the random-site model T c = 0.043 
determined in Ref. pQj. 

Further insights on the properties of the correlation 
functions can be gained by performing a scaling analysis 
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FIG. 1: Two-time charge autocorrelation functions for T — 
0.07 (a), and T = 0.03 (b). Aging effects (loss of time transla- 
tion invariance) become appreciable bellow the crossover tem- 
perature T g ss 0.05. 

of the data. We discuss separately the cases of high and 
low temperatures. 



A. T>T g 

In this temperature region the system reaches equilib- 
rium within the simulation time. We display in Fig. [51 the 
equilibrium correlation function C(t) obtained for several 
temperatures in the range 0.05 < T < 0.1. As shown 
in the inset to Fig. [2 the curves for the various tem- 
peratures collapse rather well into a single master curve 
when C(t) is represented as a function of the scaled vari- 
able t/r eq (T), where r eo (T) = exp(T /T) with T - 0.45. 
The equilibrium relaxation time thus obeys an Arrhcnius 
law above the dynamic crossover temperature. Qualita- 
tively similar results were obtained in simulations of the 
two-dimensional version of the random-site model^S^ 
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of t/t^ with /i ~ 0.7, the value for which the collapse of 
the data at large times t is the best. Repeating this pro- 
cedure for several different temperatures we obtained the 
T-dependence of the aging exponent shown in Fig. yfb). 
We find sub- aging behavior at low temperatures (/x < 1) 
in the range T <T g . When T — > T g \x decreases steeply 
to a value close to zero meaning that aging effects be- 
come negligible for T > T g ~ 0.05. The figure also shows 
the size dependence of the aging exponent. It can also 
be seen that when the linear size of the system increases 
from L = 6 to L = 10 the decay of fj, near T g ~ T c be- 
comes steeper. This makes it plausible that fj, vanishes 
(and aging stops) right at the glass transition tempera- 
ture, T c . Confirmation of this hypothesis would require a 
more detailed analysis of the L dependence of our results. 



FIG. 2: The charge autocorrelation function at equilibrium 
for T > T g . The curves correspond to T=0.10, 0.095, 0.09, 
0.085, 0.080, 0.075, 0.070, 0.065, 0.060, 0.055 and 0.050, from 
bottom to top. The inset shows the same data but represented 
as a function of the rescaled time t/r eq (T), with T eq (T) — 
exp(0.45/T). 

In previous work on the 3D model using the non-local 
dynamics described above a power-law divergence of the 
relaxation time was reported at the transition tempera- 
ture T c ^l In our case this temperature T c is of the same 
order as the dynamic crossover temperature T g for which 
our samples stay out of equilibrium during the entire sim- 
ulation time. Therefore, we have no access to the equi- 
librium critical dynamics of the model. 

B. T <T g 

This is the region in which non-equilibrium slow re- 
laxation and aging are observed. Aging effects can be 
quantified by performing a non-stationary scaling anal- 
ysis of the two-time autocorrelation functions. Experi- 
mental data in glasses are often analyzed in terms of the 
scaling formiLi^ 

C(t + t w ,t w )^F(^±^j , (5) 

where h(u) is known as the time-reparameterization func- 
tion. A commonly used form is h{u) = exp(u 1_AI /(l— /i)). 
Since this form implies an effective time scale growing 
with t w as ~ we shall analyze our data in terms of the 
simpler expression 

C(t + t w ,t w ) • (6) 

Figure |3fa) illustrates the procedure for T = 0.03. The 
inset to the figure shows C(t + t w ,t w ) as a function of 
t for three waiting times, t w = 10 3 ,10 4 and 10 5 . The 
main figure presents the same data plotted as a function 



IV. ELECTRON DIFFUSION 

The local charge correlation function discussed in the 
previous Section does not provide direct information on 
the dynamics of current fluctuations in the medium. In- 
formation on this essential aspect of the physics of the 
Coulomb glass may be obtained from an analysis of car- 
rier diffusion. 

Let Yi(t) — (xi(t),yi(t), Zi(tj) denote the position vec- 
tor of an electron at time t, where Xi, yi and Zi(t) are its 
coordinates before they are folded back into the simula- 
tion cell. We then have 

r 4 (t)=r 4 (0)+£<h-i(fc), (7) 

fc=0 

where rj(0) is the electron's initial position, Ni(t) is the 
total number of hops that it performed up to time t and 
8vi{k) is the displacement associated with the fc-th ac- 
cepted move. 

The mean squared displacement between times t w and 
t + t w is defined as 

1 M 

A(t + t w ,t w ) = — ^2{^x 2 l {t + t w ,t w )) , (8) 

i=l 

where 

A X Kt + t W ,t W )= [riit + tw \- Ti{tw)]2 , (9) 

and the angular brackets denote as before an average 
over realizations of the disorder, initial conditions and 
the thermal histories. 

Figure 01 shows the t dependence of A(t + t w ,t w ) for 
three values of the waiting time, t w = 10 3 , 10 4 and 10 5 . 
Data are displayed for two temperatures, T = 0.03 < T g , 
[FigUa)], and T = 0.06 > T g , [FigHb)]. 

It can be seen that A(t + t w ,t w ) is time translational 
invariant for the higher temperature but exhibits aging 
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T 

FIG. 3: Characterization of aging for T < T g . (a) Autocor- 
relation functions C(t + t w ,t w ) vs. the rescaled time variable 
t/t w at T = 0.030 < T g . The waiting times are t w = 10 3 
(+), tw = 10 4 (0) and t w = 10 5 (A). The inset shows the 
same data represented as a function of time, (b) Temperature 
dependence of the aging exponent fi for different system sizes 
N = L s with L = 10 (□), L = 8 (A) and L = 6 (0). 



for the lower one. Moreover, we observe that in the equi- 
librium regime T > T g the average motion is diffusive, 
A(t + t Wl t w ) = A(t) oc t, while this is noi the case in the 
aging regime T <T g . 

We discuss first the case of of high temperatures. In 
this case we can define a diffusion constant through 
A(t + t w ,t w ) — D(T)t. We attempted to fit our re- 
sults using the stretched-exponential expression D(T) ~ 
exp[-(Ti/T)' 3 ]. Using this form, a plot of T@ lnD" 1 as a 
function of T should result in a horizontal line. The inset 
to Fig. 0a) shows that this is indeed obtained for j3 ~ 1. 
The figure also shows the result of a similar plot using 
the value (3 = 0.5 expected from the Efros-Shklovski vari- 
able range hopping law. It is apparent that our data can 
not be described with this value of (3. Similar deviations 
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FIG. 4: Mean squared displacement A(t + t w ,t w ) as a func- 
tion of t for t w = 10 3 (+), tw = 10 4 (0) and t w = 10 5 (A). 
The solid line is the diffusion limit, A oc t. (a): T = 0.03 < 
T g . The dotted lines are fits of the data to Eq. 110H over the 
last decade, 10 5 < t < 10 6 . (b): T = 0.06 > T g . 

from the Efros-Shklovski law were also found in the 2D 
version of the modeli^2i for which we found (3 3/4. 27 
We now turn to the analysis of the low-temperature 
results. In this case the time dependence of A can not 
be described by a simple power law but we can still char- 
acterize the diffusion in this regime by fitting the t de- 
pendence of A(t + t w ,t w ) for our longest times (the last 
decade in t, for example) to an expression of the form 

A(t + t w ,t w ) ~f*»» T ) . (10) 

Equation l|l(J|) defines an effective diffusion exponent rj 
that depends on both the waiting time and the temper- 
ature. The asymptotic fits for T = 0.03 are represented 
by the dotted lines in Fig.^Ji where the normal diffusion 
limit, r\ = 1, is also shown for comparison. 

Figure summarizes our results for the diffusion expo- 
nent in the aging region. The temperature dependence 
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FIG. 5: (a) T dependence of the diffusion exponent rj for 
t w = 10° (+), t w = 10 1 (*), t w = 10 2 (0), t w = 10 3 (A), 
t w = 10 4 (□) and t w = 10 5 (x), from bottom to top. The 
dashed line indicates the crossover temperature T g . Inset: 
T dependence of the diffusion constant D in the equilibrium 
regime, T > T g . We plot T p log D' 1 as a function of T for 
/3 = 1 (□) and f) = 0.5 (<>)• (b) Waiting-time dependence of 
77 for T = 0.02 (+), T = 0.03 (*), T = 0.04 (0), T = 0.05 
(A), T = 0.06 (□), and T = 0.07 (x). 



of 77 for several values of t w is shown in Fig.jSJi. It can be 
seen that in the equilibrium regime, T > T g , r\ = 1 for 
all t w and T. Below T g , however, the diffusion exponent 
decreases with decreasing temperature for all values of 
t w reflecting that electron motion becomes increasingly 
sluggish. 

The dependence of r\ with t w at fixed T is shown in 
Fig. \Sjp for several temperatures. The exponent 77 in- 
creases with t w , eventually reaching the diffusion limit 
77 = 1 for long waiting times. At low temperature this 
variation is slow (approximately logarithmic) . This is yet 
another manifestation of the slow relaxation that char- 
acterizes the glassy phase of the Coulomb glass. 



The characterization of diffusion through a diffusion 
exponent is familiar in the study of random walks in ran- 
dom media were one generally finds sub-diffusive behav- 
ior (77 < 1) in those physical situations in which the dis- 
tribution of the time intervals between successive hops of 
the diffusing particle has sufficiently long tails that the 
central limit theorem no longer holds^ It would be in- 
teresting to examine the distribution of this times in the 
aging regime of our system. 



V. HETEROGENEOUS DYNAMICS 

To establish a relationship between the observed aging 
effects and the microscopic motion of electrons we ana- 
lyzed the evolution of the diffusion front. This latter is 
defined through the probability density of the squared 
displacements: 

HtiAx 2 ^,^) = -^J2\ 6 ( Ax2 - ^(t + t w ,t)) 

(11) 

where 5{x) is the usual 5- function. It is easy to show that 
for a stationary and homogeneous diffusion process the 
above distribution takes the form 



Dt 



(12) 



where $(a;) is approximately Gaussian. Hf lS thus ex- 
hibits a single peak whose position increases linearly with 
time. To compute Hi from our numerical data we must 
appropriately coarse-grain the variable Ax 2 . Since we 
found that the distribution of electron squared displace- 
ments is very broad we chose a regular coarse graining 
in log(Aa; 2 ). In Fig. EJa) we show the evolution of the 
computed Hi(Ax 2 , t, t w ) as a function of log(Aa; 2 ) for 
t w = 10 6 and T = 0.060 > T g . For this value of the tem- 
perature Hi is independent of the waiting time if t w is 
large enough. This is consistent with the observed time 
translation invariance of the charge autocorrelation func- 
tions and the mean squared displacement at this temper- 
ature. The histograms shown in the figure were scaled 
conveniently and shifted vertically by an amount log(t) 
to make their baselines coincide with the time they cor- 
respond to. 

Note that the diffusion front, located initially at Ax 2 = 
0, splits rapidly in two peaks. The position of the first 
peak is almost time independent at long times. This 
peak corresponds to squared displacements smaller than 
the average impurity distance, dQ. The center of the sec- 
ond peak increases linearly with time and its location 
coincides with the mean squared displacements at long 
times. The interpretation of these results is that the elec- 
tron dynamics is heterogeneous and characterized by the 
existence of two dynamical modes that can be clearly 
distinguished: (a) Diffusive mode: electron motion is un- 
bound and diffusive. It corresponds to metallic hopping 



7 



X 

< 




2 
log(Ax2) 



100.0 




FIG. 6: (a) Evolution of the diffusion front _ffi(Aa; 2 ), for 
t w = 10 6 and T = 0.060. Histograms are scaled and shifted 
vertically by an amount log(t) so that their baselines coincide 
with the time they correspond to. The symbols (0) represent 
the mean squared displacement A(t + t w ,t w ). The vertical 
lines are located at the positions a 2 , and L 2 , where ao is the 
mean distance between sites and L = 10 the linear size of the 
system, (b) Typical squared displacements of a selected set 
of electrons as a function of hop number. The temperature is 
T = 0.060. We show the ten first hops after t w = 10 6 steps. 
The selected graphs illustrate typical dipolar and diffusive 
electron trajectories. 



as found in Ref ,10| (b) Confined mode: electrons remain 
confined within small regions of space during the obser- 
vation time. 

Examples of trajectories of electrons of these two types 
are shown in Fig. EJb) at the same temperature and for 
the same waiting time as above. The displacements are 
represented as a function of the hop number iVj. We only 
show the ten first hops after a waiting time t w = 10 6 . 

Some of the trajectories correspond to electrons that 
hop back and forth between two sites. Ax 2 then oscil- 
lates between and &?, the distance between the sites 



involved in the motion. This fluctuating dipolar mo- 
tion contributes to most of the weight of the first peak 
in the histograms of Fig. Ela). It is important to note 
that, although this fluctuating motion appears regular 
when plotted as a function jVj, it is in fact extremely 
irregular when viewed as a function of time, since the 
time intervals between successive jumps are very widely 
distributed. The rest of the trajectories shown in Fig- 
ure HJb) are those of diffusive electrons that contribute 
to the metallic peak in the histograms of Fig. EJa). 

We turn now to the analysis of the statistics of hop- 
ping rates which is important to understand the source 
of spontaneous fluctuations in the system. To this end 
we consider the joint distribution function 

H 2 (Ax 2 ,n h ;t,t w )= (13) 

J2 t (s(Ax 2 (t + t w ,t w ) - Ax 2 )6( ni (t,t w ) - n h )^ 

where rii = iVj/ Nh is the number of hops of an electron 
i normalized by Nh = Ni/M, the average number of 
hops per electron. 

We use a regular coarse graining in log(Aa; 2 ) and rih 
to compute the corresponding histograms. These are dis- 
played in Figs. Ha) and E£b) for T = 0.070 > T g and 
T = 0.040 < T g , respectively. Both plots correspond to 
t = t w = 10 6 . 

The two dynamical modes described above can be eas- 
ily identified in the equilibrium situation [Fig. Ufa)]. The 
single peak at large Ax 2 corresponds to the diffusive 
mode while the two ridges at low values of Ax 2 corre- 
spond to the dipolar mode. Note that the hopping rates 
of electrons involved in dipolar motion have a much wider 
distribution than those of diffusive electrons. 

Fig[7Jb) shows that the two modes are still distin- 
guishable at low temperatures, when the system is out 
of equilibrium. The structure of the modes is qualita- 
tively different, however. Not only the distribution of 
hopping rates is now much broader but a small fraction 
of electrons with low hopping rate lies in between the two 
modes. The presence of these electrons, that can not be 
clearly associated with any of the modes, suggests that 
a very slow exchange of carriers between them may take 
place in the course of the relaxation. We shall further 
discuss this issue below. 

To explore in more detail the properties of electrons 
contributing to each of these modes we found it conve- 
nient to define for each electron the variable 



jr. \ _ 1 { T + t. w ,t w ) 

T=l 



(14) 



This quantity characterizes the mobility of the electron 
during a time interval of length t after t w . For a carrier 
that diffuses normally in this time span with a diffusion 
constant D, di ~ D. For an electron that remained con- 
fined in a region of linear size a during this time interval, 
di ~ a 2 hi(t)/t. Finally, for a frozen electron, i.e., one 



8 



Metallic hopping 



H 2 (Ax 2 ,n hop 




H 2 (Ax 2 ,n hop: 



FIG. 7: The joint probability distribution of squared dis- 
placement Ax 2 and normalized hopping rates H2(Ax 2 , n^) 
for T = 0.07 (a) and T = 0.04 (b). The histograms are taken 



at time t = 10 , after a waiting time t v 



10° 



that did not move at all in the time interval under con- 
sideration, dt = 0. 

We study the probability density of d defined as 



M , 

»=i » 



(15) 



As discussed above we expect f(d, t w ) to exhibit to well 
separated peaks: one at d — D(t w ), the average diffusion 
constant of the diffusing electrons at time scale t w , and 
the other at d = a(t w ) 2 \n(t)/t where a(t w ) is the typical 
size of the regions of confined motion at the same time 
scale. The width of the peaks give the dispersion of these 
quantities. 

We also introduce the cumulative distribution function 



F(d,t v 



f(x,t w )dx 



(16) 



Results are displayed in Fig. where we plot f(d,t w ), 
and its cumulative distribution function for three tem- 
peratures, T = 0.02 < T g , T = 0.035 < T g and 
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FIG. 8: Probability density f(d) (left panels) and the corre- 
sponding distribution function F(d) (right panels) for three 
temperatures, T = 0.02 (a,b), T = 0.035 (c,d) and T = 0.06 
(e,f). In all cases t = 10 6 and the waiting times are t w = 10 n 
with n = 0,1,2,3,4,5,6. The solid line in (e) indicates the 
location of the diffusion constant D(T) calculated from the 
mean squared displacements of electrons at T = 0.060. 



T = 0.060 > T g and waiting times t w = 10™ with 
n =0,1,2,3,4,5 and 6. In all cases we see the appearance 
of the two peaks referred to above. The distributions 
are stationary for T > T g for which the system equili- 
brates rapidly but they show aging in the non-equilibrium 
regime. 

A striking feature of the function / in the aging regime 
[cf. Fig. Eta,c)] is that the positions of the peaks are al- 
most independent of t w . This means that the diffusion 
constant of the "metallic" electrons is time- independent. 
The height of the peaks, however, does depend on time 
scale: as time elapses from the quench, the proportion of 
diffusing carriers diminishes while that of confined ones 
increases. This is a direct manifestation of the exchange 
mechanism that we hinted at above. The fact that the 
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location of the peaks is only weakly dependent on t w indi- 
cates that the effective mobility of the diffusive electrons 
is not much affected by the slow changes of the environ- 
ment that result from the aging process. 

The plateaus that appear in the corresponding cu- 
mulative distribution functions right after the peaks [cf. 
Figs. IHtb,d,f )] can be used to measure the relative popu- 
lations of the modes. We see a first plateau correspond- 
ing to the area of f(d,t w ) below the dipolar peak and 
a second plateau corresponding to the additional area 
below the metallic peak. It can be seen that the cumu- 
lative distribution function does not saturate to unity at 
low temperatures [cf. Figs. [Etb,d)] while it does at high 
temperatures [cf. Fig. E[f)]. The difference is due to 
the fact that, at low temperature, a fraction of the elec- 
trons remain frozen during the observation time. These 
were not counted in the numerical evaluation of the inte- 
gral in Eq. (jl fijl . Another manifestation of the presence 
of frozen carriers is the pronounced asymmetry of the 
two lower ridges in the lower panel of Fig. [7| in the zero 
hopping-rate limit. 

We can now use the height of the plateaus in 
Fig. |Sfb,d,f) to measure the populations of the differ- 
ent modes. These are represented as a function of tem- 
perature for t w — I0 6 in Fig. E^a). It is seen that the 
proportion of diffusive electrons increases with increasing 
temperature while, at the same time, that of dipolar and 
frozen ones decreases. Note that a crossover between the 
regime dominated by diffusive electrons and that domi- 
nated by confined ones is located precisely at T g . 

The waiting-time dependence of the populations is 
shown in Fig.^Jb) for two temperatures, T = 0.03 < T g 
(left panel) and T = 0.06 > T g (right panel). These pop- 
ulations are time independent at the highest temperature 
but vary logarithmically with time in the aging regime. 

It is quite tempting to try to relate these observa- 
tions to the relaxational properties of the conductiv- 
ity. Assuming that the Einstein relation holds in the 
non-equilibrium regime, the conductivity at scale t is 
a oc n f dDf(D,t), where n is the electron density and 
the integral extends over the diffusive peak. We saw ear- 
lier that the position of the peak D is time independent. 
This implies a ~ npr>D where po is the fraction of dif- 
fusing carriers. Since the latter decreases logarithmically 
with time, this simple argument predicts logarithmic re- 
laxation of the conductivity which is one of the main ex- 
perimental observations. Whether the assumptions lead- 
ing to this result are valid has to await direct computa- 
tion of the current in the aging regime, in the presence 
of an applied electric fields 



VI. CONCLUSIONS 

We have studied the relaxational properties of the 
three dimensional random-site Coulomb glass model after 
a quench from high temperature. We found a crossover 
from stationary to slow non-stationary dynamics at a 
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FIG. 9: (a) Density of diffusive (□), confined (A) and frozen 
electrons (0) during timescale t = 10 6 as a function of the 
temperature. Waiting time is t w = 10 5 . The line indicates 
the location of the crossover temperature T g . (b,c) Waiting- 
time dependence of the same densities for T = 0.03 < T g (b) 
and T = 0.06 > T g (c). The meaning of the symbols is the 
same as in (a) 



temperature T g that is very close to T c , the equilibrium 
glass transition temperature of the model. This crossover 
can be seen even in relatively small samples because of 
the exponential increase of the equilibration time with 
decreasing temperature. 

We found that at low temperature the dynamics of lo- 
cal charge fluctuations and that of current fluctuations 
show aging. In the former case, the relaxation obeys 
simple scaling laws characterized by a temperature- 
dependent aging exponent n(T). Analysis of the tem- 
perature and system-size dependence of /^(T) suggests 
that in the thermodynamic limit the observed crossover 
at T g ~ T c becomes a real dynamic transition that occurs 
precisely at T c . 

The analysis of the properties of diffusion fronts re- 
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vealed that the dynamics of carriers is heterogeneous as 
it was found in other glassy systems i2i We found that, for 
each timescale two classes of electrons may be identified, 
those that have diffusive motion during the observation 
time and those whose motion in the same time interval 
remains confined. Only electrons belonging to the former 
class contribute to the dc conductivity while the others 
only contribute to the dielectric screening. 

In the region of low temperatures where aging is ob- 
served electrons are slowly exchanged between these 
two modes with the consequence that the population 
of metallic electrons decreases logarithmically with time 
without appreciable change of the their diffusion con- 
stant. This provides a plausible explanation for the loga- 
rithmic relaxation of the conductance after a quench that 
was observed experimentally. 

We believe that the local microscopic dynamics used 
here, which is a realistic description of hopping processes 
in experimental systems, plays an important role in the 



phenomena that we observed. This type of dynamics fa- 
vors the appearance of local effective constraints that re- 
late our model to kinetically constrained models in which 
slow dynamics arises from restrictions on the allowed 
transitions between configurations. 
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